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I. INTRODUCTION. 

Topological insulators represent a nev^ state of mat- 
ter where the topology stabilizes some extremely robust 
properties such as the emergence of conducting metallic 
edge or surface states in spite of the presence of im- 
purities or defects. Topological insulators have been 
theoretically predicted in 2- and 3-dimensions and now^ 
v^e have several concrete materials that sYvow some of 
the predicted properties. 12^1 For the 3-dimensional topo- 
logical materials, the surface states have been mapped 
using angle resolved photoemission spectroscopy and 
they show an odd number of Dirac cones, in perfect 
agreement with the theory.SEI Fiowever, the transport 
measurements have revealed that, so far, all the samples 
display a metallic character in the bulk. As such, the 
transport experiments have taken a central stage in the 
research on topological insulators. 

The literature abounds with high quality exp erimen- 
tal transport data for topological insulators.'^'^^^ There 
are detailed reports about the behavior of the transport 
coefficients with the temperature. Magneto-electric mea- 
surements are also available, with maps of the conduc- 
tivity tensor as function of magnetic field strength and 
electron density, which is controlled by gate voltages or 
by doping. There are also studies done at finite frequen- 
cies and maps of the transport coefficients as functions 
of films' thickness and disorder strength have been re- 
ported. The materials have been progressively tuned, to 
a point where the contributions to the transport from the 
bulk and surface are comparable. The available experi- 
mental data could be used as a window into the micro- 
scopic properties of these materials, if efficient quantita- 
tive theoretical analyses could be developed and applied 
to the real materials. Such analyses will have to include 
the disorder and the magnetic fields, a task that at first 
sight may seem extremely difficult. 

We will argue here that actually that is not the case: 
The disordered systems under magnetic field can be an- 
alyzed very much like we analyze the translational in- 



variant systems. Let us take a few lines to explain what 
we mean by this. For translationally invariant systems, 
the response and the correlation functions can be con- 
veniently computed in the dual fc-space. The fc-space 
analysis consists of two parts: 

PI: Derivation of closed-form expressions. 

P2: Numerical evaluation of the formulas. 

For translational invariant systems, no matter how 
complex the correlation function, one can easily derive 
closed-form expressions, which typically involve inte- 
grations and derivations of ordinary functions defined 
over the Brillouin torus. Even though these formulas 
are formal, in the sense that one still needs a computer 
to evaluate them, they pretty much fulfill our idea of an 
analytic solution because the expressions are transpar- 
ent enough to enable qualitative understanding and, on 
the quantitative side, the formulas can be numerically 
evaluated without much effort. 

The analog of the fc-space calculus for aperiodic sys- 
tems is the noncommutative Brillouin torus and its non- 
commutative calculus developed by Bellissard and his 
collaborators.^^ This formalism has been used already 
in the literature to derive closed-form expressions for 
the Kubo formula of transport^^ electric polariza- 
tion and orbital magnetization,^! all in the presence of 
disorder and magnetic fields. We have used the non- 
commutative calculus to compute topological invariants 
for disordered topological insulators .^^ ^"^ Fiastings and 
Loring have also used the noncommutative setting to 
derive and to efficiently compute new invariants for dis- 
ordered topological insulators-S^'I^These applications of 
the noncommutative calculus pretty much demonstrate 
that any response or correlation function, that can be 
written in closed-form in the fc-space, can be translated 
into a closed-form noncommutative expression which 
incorporates the effect of disorder and magnetic fields. 
As such, the first part (PI) of the analysis is already in 
place. In a previous work,^^ we have made progress on 
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the second part (P2) of the analysis, namely on how to 
quantitatively evaluate the non-commutative formulas. 

For translational invariant systems, the fc-derivatives 
are evaluated using finely tuned finite-difference algo- 
rithms and the integrals in the fc-space are evaluated 
using Riemann sums, both done for various samplings 
of the Brillouin torus. This leads to approximate results 
that converge to the exact result as the sampling of the 
Brillouin torus becomes finer and finer. If the integrands 
in these formulas are analytic functions of fc, as it is the 
case at finite temperatures, the convergence happens ex- 
ponentially fast. This is an important characteristic of 
the fc-space calculus because it enables extremely accu- 
rate calculations even with modest computational ef- 
forts. Note that in a first-principle calculation, one often 
deals with hundreds of energy bands so even for fc-space 
calculus the efficiency is a big issue. 

Now, when disorder and magnetic fields are present, 
we have to switch to the noncommutative Brillouin torus 
and its noncommutative calculus, which are defined in 
the strict thermodynamic limit. These structures don't 
have immediate analogs at finite volumes, but as it was 
elaborated in Ref. |5T} one can define an approximate, 
finite-dimensional non-commutative Brillouin torus and 
an approximate non-commutative calculus. These ap- 
proximate structures can be implemented and manipu- 
lated on a computer and, as such, the non-commutative 
formulas can be numerically evaluated, though in an 
approximate way. One key result of Ref. 51 is that the 
approximate results obtained in this way converge expo- 
nentially fast to the exact ones as the dimensionality of 
the non-commutative Brillouin torus is increased. Note 
that an inverse power law convergence is not sufficient 
for exploring the typical interesting questions arising for 
disordered systems. Hence, in some sense, Ref.'Sl gave 
a solution for the second part (P2) of the analysis. 

The structure of the present paper is as follows. 
In the first part we give an introduction of the non- 
commutative Brillouin torus and its non-commutative 
calculus, and present a formal derivation of the non- 
commutative Kubo formula, closely following Ref. '40] 
Even though this material has been reviewed with other 
occasions (see for example Ref. 52), an intuitive exposi- 
tion using a language that is a little more familiar to the 
condensed matter theorists could be of interest. Further- 
more, we thought that the readership will welcome a 
format where the non-commutative framework, the nu- 
merical algorithm and the applications are all presented 
in one place. 

In the second part we present the approximate finite 
volume non-commutative Brillouin torus and its non- 
commutative calculus, together with the emerging ap- 
proximate Kubo formula. The formula is then broken 
down to an explicit expression which can be straightfor- 
wardly implemented on a computer. 

In the third part we present an application to a model 
of a disordered 2-dimensional quantum spin-Hall insu- 
lator without edges. We map the bulk conductivity ten- 



sor as function of Fermi energy (Ef ), disorder strength 
(W) and temperature. We report several convergence 
tests. Based on this calculations, we identify the metallic 
phase and we map the phase diagram of the systems 
in the (£f , W) plane. The results give a direct confir- 
mation, via the computation of the conductivity tensor, 
that strong disorder closes the insulating gap and drives 
the system into a metallic phase, and then into a topo- 
logically trivial phase. When we place the model in 
the trivial phase, we find that the disorder alone can 
drive the system into a metallic phase. The phase di- 
agrams are found to be in very good agreement with 
the ones reported in previous studies. Furthermore, the 
phase diagram derived from the transport simulations is 
compared with the phase diagram derived from a level 
statistics analysis and good agreement is found. 

In a second application, we turn the magnetic field 
on, and we point out a markedly different behavior of 
the resulting Hofstadter spectra for topological versus 
non-topological systems. We then map the resistivity 
tensor as function of Fermi level at fixed magnetic field. 
We demonstrate that the algorithm is able to resolve the 
expected Hall plateaus. 



II. DISORDERED LATTICE MODELS IN THE 
PRESENCE OF MAGNETIC FIELDS 

We will restrict the discussion to 2-dimensional lat- 
tice models without edges. This setting is not restrictive 
because it covers the bulk 2-dimensional Quantum spin- 
Hall insulators and the films of 3-dimensional insulators 
(all transport experiments are carried on films). As we 
shall see, there are many interesting questions arising 
even for bulk 2-dimensional Quantum spin-Hall insula- 
tors, which are investigated in this study. 

The Hilbert space of a 2-dimensional lattice model is 
spanned by functions ip defined over the 2-dimensional 
lattice with values in C^, where D is the number of 
quantum states per site. In the clean limit, D is also equal 
to the number of bands of the model. We will consider a 
general lattice Hamiltonian with on-site disorder: 

(H,t/;)(n) = E e-'^^-tm-nil^{m) + Wcbnij^in). n) 

m 

The phase factor 

^-i(pnm - ^-incpinAm) ^2) 

encodes the effect of the magnetic field via the Peierls 
substitution,^^ where cp is the magnetic flux per repeat- 
ing cell, measured in the units of flux quantum (po = hje, 
and n Am = nim2 - niin2. The D xD matrix encodes 
the hopping amplitudes from a site to its neighboring 
site situated at the relative distance fc, and a)„ is a di- 
agonal D X D matrix with independent random entries 

uniformly distributed in the interval -|/ ^ • 



3 



The collection of the matrix amplitudes {cbn}nez? can 
be viewed as a point co in the space Q defined as: 



The space Q will be equipped with the probability mea- 
sure: 

dco = U„,z^Ua=^dco"„. (4) 
The disorder average is then given by the integral 

There is a natural action of the additive discrete 
group on Q given by: 

: Q ^ Q, (tmco)^ = a)n+nf (5) 

The measure dco is invariant and ergodic relative to the 
action of this group. 

We would like to point out that accurate lattice models 
exist for most of the topological materials. They can be 
derived empirically or from first-principle calculations 
by following, for example, the methods presented in 
Ref . 54 These models are believed to be well suited for 
the transport simulations. 

III. THE NONCOMMUTATIVE FRAMEWORK 

For periodic solids, the correlation functions can be 
cast as closed formulas involving the classic integro- 
differential calculus over the Brillouin torus. These for- 
mulas give the desired answers directly in the thermo- 
dynamic limit. Furthermore, by examining the degree 
of smoothness for the functions entering these explicit 
expressions, one can easily understand if the formulas 
are well behaved and how to evaluate them numerically, 
i.e. how to proceed with the discretization of the classi- 
cal Brillouin torus and what finite-difference scheme to 
use for the fc-derivatives. Similarly, the noncommutative 
calculus provides an extremely convenient framework to 
carry the calculations for aperiodic solids and ultimately 
to derive closed formulas directly in the thermodynamic 
limit. Furthermore, like in the periodic case, the non- 
commutative calculus enables one to understand when 
these formulas are well behaved and how to compute 
them numerically. 

A. The Abstract Algebra of Observables 

For periodic solids, the integro-differential calculus 
goes over the algebra of functions defined on the clas- 
sical Brillouin torus and let us recall a classic result in 
set topology which says that a topological space can be 
reconstructed from the commutative algebra of contin- 
uous functions defined over that space. As such, the 
geometric Brilloin torus and the calculus defined over it 



can be defined in a pure algebraic setting. This is exactly 
the kind of shift of reasoning that is needed when deal- 
ing with aperiodic solids because, when disorder or a 
magnetic field with irrational flux per repeating cell are 
present, the geometric Brillouin torus loses its meaning 
but the commutative algebra of functions defined over 
the classic Brillouin torus can be easily adapted to the 
situation. This algebra is replaced by a noncommutative 
C'-algebra and, in the spirit of what was said above, the 
resulting noncommutative C -algebra can be rightfully 
called the noncommutative Brillouin torus. 1^ 

We now describe this algebra. Its elements are func- 
tions defined on Q x with values in the set Mdxd of 
D X D complex matrices. We will use lower capital letters 
like f, gorh to refer to the elements of the algebra. The 
addition rule for the algebra is: 

(/ + ^) = /(^/ ^) + g{(^' (6) 

where the last " + " operation is the ordinary matrix 
addition. The multiplication rule is: 

(/ * g){co, n)=Y^ f{co, m)g{t-Jco, n - (7) 

where the product between / and g appearing on the 
right hand side is the usual matrix multiplication. The 
algebra has a unit element defined by: n) = 6n,o- 

The link between the abstract algebra defined above 
and the operators acting on the quantum states is simple. 
Each element / generates a covariant family of operators 
TZojfr which act on a quantum state via the following 
formula: 

{{nM) (") = £ fi^n'^' « - (8) 

For example, the Hamiltonian in Eq. [T] is generated by 
the element: 

h{a),n) = tn + W5n^oCOnr (9) 

as one can easily verify that: Ucoh = H^. 

The map defines a representation of the algebra in 
the space of operators acting on the quantum states: 

nUf-g) = M){nc.gy (10) 
The covariant property means: 

Ua{n^f)U~a' = ^Ucof, (11) 

where Ua are the magnetic translations on ^^(Z^): 

(Uail^) (n) = - a). (12) 

The algebra defined above can be endowed with the 
structure of a C'-algebra, a fact that allows one to define 
a fine spectral theory and a functional calculus, that is, a 
natural framework to define functions of an element. For 
this we need a norm with the special property 11/ g|| < 
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ll/llllgll and an involution (usually called a '^-operation) 
f^f such that II = ll/lp. If II • ir denotes the usual 
norm for operators acting on the quantum states, then 
the following norm and the star operation: 



= sup ||7T^/ir, f{a), n) = /(t_„a;, -n)^ 



(13) 



have those properties. If we consider only those ele- 
ments for which the norm in Eq. 13 is finite, then the 
algebra of observables becomes a c^-algebra which is 
denoted as It is this algebra that is called the non- 
commutative Brillouin torus. 



B. Spectrum, the resolvent function and the analytic 
functional calculus 



The eigenvalue problem: 

Hip A = Ai/^A (14) 

and the functional calculus 

^{H) = La^{A)\iP^){iP^\ (15) 

for ordinary self -adjoint matrices can be formulated in 
a purely algebraic language. In fact, both concepts can 
be naturally developed and generalized in the abstract 
setting of C -algebras, without making any reference to 
the eigenvectors. A good reference for this^topic is the 
short course in spectral theory by Arveson.l^ 

Let us first discuss the notion of the spectrum of an 
element / from the algebra which generalizes the set 
of eigenvalues. The points z of the complex plane for 
which z - / is invertible (in ^) form an open set, called 
the resolvent set of /. The usual notation for this set is 
p{f). The spectrum of the element / is the complement of 
the resolvent set: C - p(/), and this set is usually denoted 
by cr(/). For any / in a C-algebra, the spectrum a{f) is a 
non-empty compact subset of the complex plane. There 
are two particular classes of elements that we want to 
mention: 1) the self-adjoint elements: 



r=f> 



(16) 



whose spectra are confined on the real axis, and 2) the 
unitary elements: 



(17) 



whose spectra are confined on the unit circle. For exam- 
ple, the Hamiltonian h defined in Eq. [9]is a self -adjoint 
element while the time evolution generated h is a one- 
parameter group of unitary elements. 

The inverse of z - / will be denoted by (z - /)~^ . It is an 
element of Ji. When viewed as a function of z, (z - f)~^ 
is called the resolvent function of /. It is an analytic 
function of z on the resolvent set p{f) with values in 
Given a function 0(z), analytic in a neighborhood of the 



spectrum of /, one can define the function O of / through 
the formula: 



(18) 



where C is a contour confined to the analytic domain of 
0(z) and circling the spectrum of /. The above formula 
provides a functional calculus, i.e. a morphism between 
the algebra of analytic functions in a neighborhood of 
a{f) and the algebra Ji, that is: 



((i,.(D')(/) = o(/)*om 



(19) 



for any O and two such functions. As such, any 
algebraic identity for ordinary functions, belonging to 
the class allowed above, translates automatically into an 
identity for the same functions but with z replaced by /. 



C. The noncommutative integro-differential calculus 

The commutative algebra of ordinary functions de- 
fined over the classical Brillouin torus has been replaced 
by the non-commutative C-algebra The correlation 
functions for translational invariant systems are com- 
puted via integrals of the form: 



/ dk tr{f (fc)}. 



(20) 



where tr is the trace over Mdxd space. The operation 
in Eq. |20l can be seen as a linear functional defined on 



the space of functions defined over the classical Brillouin 
torus with values in Mdxd- This linear functional has 
two special properties. It is cyclic: 



/ dk tr{F{k)G{k)] = J dk tr{G{k)F{k)] 
and it is positive: 

/ dk tr{f (/c)+f (fc)} > 0. 



(21) 



(22) 



These two properties makes the linear functional of 
Eq. 20 into a generalized trace. 

Now on it is still possible to define a bounded 
linear functional that replaces the fc-integration. This 
functional is: 



nf) = f^dco tr{/(a;,0)}. 



(23) 



where, like before, tr is the trace over Mdxd space. The 
linear functional defined in Eq.|23]is cyclic and positive: 



T{f*g) = ng*f), T{f*n>o. 



(24) 



This makes T into a generalized trace and T replaces the 
fc-integration. The existence of a generalized trace is also 
extremely useful because the algebra Ji can be endowed 
with a scalar product: 



(f,g) = nrg) 



(25) 
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which makes ^ into a Hilbert space. As we shall see, 
this is important because the linear maps acting on the 
space of operators then become linear maps on a Hilbert 
space and we do know how to manipulate linear maps 
on Hilbert spaces. 

IS 



The noncommutative integration defined in Eq. 23 



natural in the sense that in the absence of disorder it re- 
duces to the ordinary fc-integration. Furthermore, like 
the classic integration over the Brillouin torus, the non- 
commutative integration defined in Eq.[23]gives the cor- 
relation functions directly in the thermodynamic limit. 
More precisely, if f , G, . . ., are the operators acting on the 
quantum states generated by /, g, • • ., then: 



Vol{A) 



TrA{fG...} = r(/>^g^...), 



(26) 



where the trace at the left hand side is taken only over 
the quantum states inside the box A. The identity of 
Eq.l^not only gives a close formula for the correlation 
functions for disordered systems under magnetic fields, 
but also provides a convenient way to determine when 
a correlation function is well defined in the thermody- 
namic limit. We will have an entire discussion of this 
aspect for the conductivity tensor. 

For translationally invariant systems, most if not all 
the correlation functions of interest are of the form: 



/dfctr|5,«'Oi(H(fc))<9^^02(H(fc)). 



(27) 



where ^i{H{k)) are some functions of the Hamiltonian 
and (9^' are the fc-derivations (possibly of higher power) 
on such functions. Thus, the fc-derivations are important. 
On the algebra one can define a set of automorphisms 
to replace the fc-derivations: 



dkjF{k) ^ {djf){a),n) = injf{a),n), 



(28) 



where ny is the ;-th component of n. The derivations 
defined above are natural in the sense that in the absence 
of disorder, they reduce to the ordinary fc-derivations. 
Some of the classic rules in integro-differential calculus 
still apply, such as didj = djdi, di{f^g) = {dif)^g + f^{dig) 
(Leibniz rule) or T{0{h)diO\h)) = 0. 

We end by pointing out that the representation of the 
derivation on the quantum states is: 



(29) 



As such, the element of ^ that generates the charge- 
current operator is: 



y = -v/z. 



smce: 



where e = -1 is the electron charge. 



(30) 



(31) 



IV. THE NONCOMMUTATIVE KUBO FORMULA 

We now have all the rules of calculus and we can pro- 
ceed with the derivation of the noncommutative Kubo 
formula. We closely follow Refs. 39 and 40 and we will 
present only the main steps of the derivation. The inter- 
ested reader can consult the ample review article from 



A. The coherent time evolution 

The Hamiltonian h itself defines a derivation, i.e. a 
linear map on the algebra ^ that satisfies the Leibniz 
rule,^^ through the Liouvillian: 



£h[f] = i{h*f-f*h). 



The equation: 



dtu(t) = -£h[u{t)l m(0) = 1, 



(32) 



(33) 



defines a one parameter unitary flow u{t) over 
which implements the time evolution in the Heisenberg 
picture. Explicitly: 



-ith 



(34) 



In the presence of a uniform electric field E, the Hamil- 
tonian becomes (here we sete = -1): 

hE=h + E'V (35) 

and the quantum time evolution is generated as: 

UE{t)f = e-'^^Ef^ (36) 

with £h^=£- EV. 

B. Scattering processes and the dissipative time evolution 

Decoherence is introduced by random scattering 
events, such as the electron-phonon scattering. The scat- 
tering matrices for various dissipation mechanisms are 
known explicitly.^^ Here, however, we will proceed at 
a formal level and assume that the scattering matrix is 
known and given. Later on, we will adopt the relaxation 
time approximation. 

If the scattering matrix for a scattering event at time t 
is denoted by Wt^ the instantaneous time evolution from 
initial time t = Oto some arbitrary time t takes the form: 

UwA^) = UE{t - tn)Wt^UE{tn " tn-l)Wt^_, • • • Wt.UEih). 

(37) 

The time sequence {fz}z=i,2,... is assumed to be generated 
by a Poisson random process with an average frequency 
1/t. The scattering matrices w can fluctuate from a scat- 
tering event to another. The average of We(0 over the 
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collision times and the fluctuations of w's, denoted by 
UE{t), can be computed explicitly, and the result is:^^^° 



UE{t) = e 



(38) 



where F = (1 - zZ))/t with w being the average of w over 
its fluctuations. 

The evolution z^e(0 includes the dissipation and is no 
longer unitary However, in the absence of an electric 
field, the dissipative time evolution must leave the ther- 
mal equilibrium state un-altered, and this requires that 
r and X/z commute. 



C. The conductivity tensor 

We assume that the electric field was turned on at f = 0, 
when the system was still in its thermal equilibrium state 
described by the density matrix po = ^fd(^)- Here, Ofo 
is the Fermi-Dirac distribution function corresponding 
to a temperature T and Fermi level Ef. After the field 
was turned on, the density matrix evolves according to: 
p(0 - Uw,E{t)po- As such, the instantaneous expected 
value of the charge current density is: 

/z.,E(0 = T(y^u^,E(Opo). (39) 
The average over the collision processes leads to: 

jE{t) = T{j*UE{t)po), (40) 

and the time average: 

(7E) = liml Jdt'M') (41) 
can be computed as it follows: 

(7E) = lim6 f^^dte-'^jEit). (42) 



With the explicit expression of /^(f) from Eq. 40 and the 
explicit formula for z^e(0 from Eq. 38 the integral over t 
can be computed explicitly leading to: 



(Je) = lims^o 6 r(y - (6 + F + £h,r'po) • 



(43) 



Since there is no current in the absence of the electric 
field, we must have: 



(44) 



and we can subtract such null term from Eq. 43 which 
then becomes: 



(jE) = lim6T(j-{6 + Y + £h,r 
o{EV)o{6 + T + £f,)-^po). 



(45) 



Since JlhPo = and Fpo = (recall that X/z and F com- 
mute), the action of the map seen in the second line of 



the above equation reduces to just 6 ^po. Hence, we can 
take the limit 6 ^ to finally obtain: 



(7e> = -T({Wh) * (r + X,,)-i(EVpo)) . 



(46) 



The non-linear conductivity tensor can be easily read 
from above: 

aij(E) = -T {{dih) * (T + £h, )-^^/Ofd(/j)) • (47) 



Its limit as E — > gives the linear conductivity tensor: 



Oij = -T({dih) * (r + £hr'di^m{h)), 



(48) 



and this is the non-commutative Kubo formula. 

In the relaxation time approximation, which we adopt 
from here on, the map F is replace by the identity map 
times a coefficient l/Trei- The relaxation time Trei is an 
empirical parameter, which is assumed as given. The 
Kubo formula in the relaxation time approximation be- 



comes: 



Oij = -T({dih) ^ (1/Trel + X/z)-'^;Ofd(/z)) • 



(49) 



This is the expression that will be used in our present 
transport simulations. For this simplified expression, 
one can easily show that all the entries are well be- 
haved. In particular, the inverse (l/xrei + X/z)~^ exists 
in since the map z'X/i is a self-adjoint operator in the 
Hilbert space defined by the scalar product of Eq. 25 
hence the operator X/z does not have spectrum at -l/xrei 
(or in other words, -l/xrei belongs to its resolvent set 
hence l/xrei + -Lh is invertible). Also, the derivation of 
Fermi-Dirac function of h belongs to the algebra J[ due 
to its rapid decay to infinity. As such, all the entries 
in Eq. 49 are well defined hence the non-commutative 
Kubo formula takes finite values and is numerically sta- 
ble. We should point out that without the noncommu- 
tative formalism, the most one could do was to express 
the conductivity tensor as a thermodynamic limit (like 
in Eq.|26l whose existence and stability would have been 
very dmicult to establish. 



V. AN OPTIMAL KUBO FORMULA AT FINITE 
VOLUMES 

The exercise from the previous section, we hope, is 
fairly convincing in showing that the linear response co- 
efficients, in general, can be expressed as compact and 
transparent formulas. In this section we discuss how 
to efficiently evaluate such formulas, which are written 
directly in the thermodynamic limit and the thermody- 
namic limit is not accessible to a computer. Once we com- 
plete this numerical aspect, we will have a rigorous and 
practical formalism to investigate the linear response co- 
efficients of disordered systems under magnetic fields. 
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A. The C-algebra over a torus 

We consider a finite square lattice which is wrapped 
into a torus. The torus is generated by two discrete circles 
T = X S^. We specifically require that each discrete 
circle contain 27V + 1 points. On this torus we pick an 
arbitrary point and call it the origin o. We also introduce 
a coordinate system on T by naturally un-rolling the 
torus onto the points {-N, . . . , A/'}^ c Z^, such that the 
coordinate of the origin is 0. The coordinates of a point 
p eT will be denoted by tip. The total number of nodes, 
equal to (27V + l)^ will be denoted by |T|. 

We now define the group of rotations, which replaces 
the group of translations on Z^. Given a point p of the 
discrete torus, we can imagine a succession of rigid rota- 
tions of the 5^ circles, that rotate the torus until the ori- 
gin o reaches the position where the point p was located 
before the rotations. All the other points rotate rigidly 
with the origin, hence this action defines a map Xp on the 
torus. It is clear that Xp{o) = p, and that Xp is independent 
of the specific sequence (which is not unique) of rota- 
tions used to bring the origin at position p. The rotations 
satisfy the following commutative group relations: 

^oXp = Xr^p =XpOX^= Xr^^. (50) 

We now define the equivalent Q space for the torus, 
which is denoted by Qaa- Let d) = {a)p}peT be a sequence 
of diagonal D xD matrices with identical independent 
random entries uniformly distribution in [-1/2,1/2]. 
Then d) can be viewed as a point of the probability space 
= [-1/2,1/2]^™, which will be endowed with the 
probability measure: 

dcb=UpeTna=ldcOr ^^^^ 

The rotations r induce a group of automorphisms on Qyv: 

XtjO) = {cdx^p}peTr (52) 

whose actions leave the measure dd) invariant. 

The C-algebra Ji^ over the torus is defined as follows. 
The elements are functions / defined on Qyv x T and 
taking values in Mdxd- The law of composition is: 

(/^g)(^^/P) 

= E /(a),^7)g(r-id;,r-V)e^'^^("^^"'?). ^^^^ 

We need to define a norm and a ^-operation with the 
proper characteristics mentioned when we discussed the 
C-algebra As before, the norm will be introduced via 
the operator representations. Each element / induces 
an operator acting on ^^(T) x C^, the square summable 
sequences defined over the torus with values in C^: 

(ificoDcp) ip) = L /(V^a;,r^i^)e-^KA-.)^(^). ^^^^ 



The map fto) provides a representation of the algebra JHt, 
i.e. 

^df"" g) = {Ticbf){na,g), (55) 

if and only if the flux (p (in the units of (po) takes a quan- 
tized value: 

2 

^ " 27V + 1 ^ 

which, from now on, it will always be assumed to be 
true. We mention that the quantization condition for the 
flux is in line with Zak's finding that the magnetic trans- 
lations accept finite representations only if the above 
quantization is satisfied. 

In this conditions, one can define the norm: 

ll/ll = sup WfiafW (57) 

where this time the primed norm denotes the operator 
norm on ^^(T) X C^. This norm has the desired property 
that 11/ gll < ll/ll llgll- Furthermore, a ^-operation can be 
defined: 

r{(b,p) = f\x-^cb,x-^o)\ (58) 

so that the n representation and the ^{-operation satisfy 
the essential relation: 

fto-f = (ft<5/y. (59) 

As such, the norm defined in Eq.|57|has the fundamental 
property: 

llAf II = 11/11', (60) 
which makes into a C^-algebra. 



B. An approximate integro-differential calculus over the 
torus 

We now introduce the generalized trace (the integra- 
tion) over It is defined by: 

'^t(/) = i^^d(bir[f{(b,o)]. (61) 

This linear functional satisfies all the required properties 
of a generalized trace, i.e. positivity and cyclicality. The 
trace can be computed via the equivalent formula: 

TT{f)=^Ja^d&lr\ns,f]. (62) 

The differential calculus defined over ^ will not work 
over the torus because ifipfip) does not close continu- 
ously at the boundaries of the coordinate system. What 
we can do is to define an approximate differential cal- 
culus, and that is done as follows. Let x : Sl^ ^ ]R be 
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a continuous function such that (Up is the coordinate of 
the point p eS^^): 



\x{p) - Mpl = if \np\ < N/2 



(63) 



and: 



Wp)|<|np|forallpG5^. (64) 
Then the formula for the approximate derivations is: 



0if){p) = Hpdfip). 



(65) 



This formula acts like the exact derivation on functions 
which takes non-zero values only around the origin 
o. The elements entering in the Kubo formula, like 
djOpoWr are concentrated near the origin and have a 
fast exponential decay away from the origin. So the er- 
rors introduced by the approximate derivations decay 
exponentially fast with the size of the torus. 

The operator representation of the derivation is as fol- 
lows. Let O be the discrete unit circle in the complex 
plane defined by the solutions of z^^^^ = 1 and let: 



(66) 



be the discrete Fourier decomposition of the function x. 
Then: 



(67) 



where Xi is the z-th coordinate operator: {xi(p){p) = 
np.(p{p) for (p a quantum state over the torus (i.e. in 
^2(T) X C^). 

In the numerical simulations, the function x{p) will be 
chosen in the same way as in Ref . 4J, where we computed 
the Chern number using the non-commutative calculus. 



C. The approximate Kubo formula on the torus 

If /z is a short-range Hamiltonian in then it also 
defines a Hamiltonian h in Ji^: 



h{d),p) = h{6j,np), 



(68) 



provided the torus is large enough. Note that d) e Oj^ 
is also part of the space Q. We now can write down the 
approximate Kubo formula on the torus: 



dij = -Tt ({dih) ^ (1/Trel + £j,)-'djOMh)) . 



(69) 



The above expression is self-averaging, so the fluctua- 
tions with d) rapidly disappear as the size of the torus is 
increased. The most important fact about the above for- 
mula is the following error estimate that was established 
inRef.ET] 



\dij - Gijl < Qe" 



which tells that Eq. 69 converges exponentially fast 
to its thermodynamiclimit. Above, ^ is any con- 
stant between the bounds < ^ < sinh~^(K/2d/z), with 
K = min{l/2Trei, TifcT} and h being the largest hopping 
amplitude of h. The constant is fully identifiable and 
increases with ^. 

The formula in Eq. |69] can be directly evaluated us- 
ing the standard linear algebra routines, but replacing 
the matrix operations with the composition rules and 
the norms corresponding to algebra Ji^, the latter be- 
ing needed to evaluate the remainders and to advance 
the iterations (hence the norms introduced above are 
not just for formalities). In the present work, however, 
we adopted a more traditional way even though it is 
probably not very efficient. It is based on the following 
observations. If l€a,(pa}i^Y^\ is the eigensystem for the 
Hamiltonian fc.-.h: 



(71) 



then Eq. 69 can be cast in the following equivalent form: 

D|T| 



((/),|7I^((9/z)|(/)^)((/)|,|7I^((90FD(/^)|(/)a) 

1/t + i{ea - €b) 



(72) 



This expression may look very similar to the usual Kubo 
formula at finite volume, already used in the previous 
simulations for disordered systems,!^^ but there are a 
few major differences. Since the derivation d has a more 
involved expression in our case, the derivation of OFD(/i) 
cannot be processed any further. As such, our expression 
doesn't look like a current-current correlation function. 
We can bring our formula to such form if we use a low 
approximations for d but that will destroy the exponen- 
tial convergence. 

At the end of our theoretical exposition, we want to 
stress again that the approximate formula Eq. [69] and 
the error estimate Eq.[70| should be considered together, 
because one of them gives a practical way to do computa- 
tions and the other really tells what is being computed. 
The C -algebraic framework has been instrumental for 
obtaining both results in Ref . |5l] We have also insisted 
on presenting the C -formalism because then one could 
follow a similar philosophy to derive and compute other 
important linear response coefficients, such as the spin 
transport coefficients. 



VI. TRANSPORT SIMULATIONS FOR A DISORDERED 
QUANTUM SPIN-HALL INSULATOR 

A. The model 



We will work with the Bernevig-Hughes-Zhang 
model,^ which fits the 'Tow energy" band theory of the 
(70) clean HgTe/CdTe wells. In the fc-space, this effective 
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Hamiltonian takes the form: 



Kk) m 



(73) 



where h{k) = e{k) + d{k) • cr, with o = (cTj, ct^, ct^) encod- 
ing the Pauli's matrices, and r(A;) is a S^-nonconserving 
interaction. The minimal form of the r(fc) matrix is:^^ 



r(fc) = zA| 



sin kx - i sin ky 

sin kx + / sin ky 



(74) 



The experimentally measured energy bands, in the prox- 
imity of the F-point, can be captured by the following 
expression for d{k): 

d = {As\nkx,Asinky,M - 2B{2 - coskx - cos ky)). (75) 

The Hamiltonian Hq displays a topological phase if 
0<M/B<4 and 4<M/B<8 with the insulating gap clos- 
ing at M/B=0, 4 and 8, and a topologically trivial phase 
if M/B<0 or M/B>S. In our simulations, the parameters 
will be fixed at: A = 1, B = 1, C = 0, D = 0, A = 0.5, and 
€{k) will be also set to zero. The phase diagram of the 
model, with exactly these same parameter values, has 
bee investigated in Refs. 62 and 45 We will compare the 
phase diagrams derived from the transport simulations 
with the phase diagrams from these two references. 

The model can be realized in real space using a square 
lattice with four quantum statesper site. Explicitly, the 
matrices ?„, to be plugged in Eq.n] are given by: 



^(1,0) = ^, 



(-1,0) 



§ai + B(J3 

_ A 
2 



A 
2 



^(0,1) - ^, 



(0,-1) 



|(J2 + Bos 



^^3 



-|(J2 + B(J3 



and 



_ (M-4B)(J3 
^(O'O) "I (M - 4B)(J3 



(76) 



(77) 



(78) 



The onsite disorder is introduced as discussed in the 
second section, with the restriction that it needs to pre- 
serve the time reversal symmetry. For this, the ran- 
dom amplitudes that enter on the diagonal of the ma- 
trix cbn (the rest of the entries are zero) must satisfy: 

id)n)n = {^n)33 and (cOn)!! = i^n)U' 



B. Simple test results 

We show here the most straightforward test calcula- 
tion, with the disorder and magnetic field turned off. 
In this case, the conductivity is given by the traditional 
Kubo formula which can be evaluated in the fc-space 
where we can increase the sampling of the Brillouin 





Gil (50x50 lattice) 


Gil (Exact) 


Error 


0.000 
-0.368 


0.013 


0.013 
0.013 





0.013 


0% 




1 0.013 


0.013 


0% 


-1.105 


166.284 


162.006 


2% 




[ 390.475 


389.231 


0.3% 


-1.842 


428.067 


429.085 


0.2% 


-2.211 


444.573 


443.920 


0.1% 


-2.579 


460.576 


459.891 


0.1% 




1 448.714 


448.550 


0.1% 


-3.316 


421.281 


422.249 


0.2% 


■■ 


P 387.195 


385.555 


0.4% 


-4.053 


341.874 


340.844 


0.3% 


-4.421 


292.086 


289.548 


0.8% 


-4.790 


232.352 


232.622 


0.1% 


-5.158 


171.330 


170.758 


0.3% 


-5.526 


104.282 


104.477 


0.2% 


-5.895 


34.388 


34.187 


0.6% 



TABLE I. The diagonal conductivity Gh computed via Eq. 72 
for different Fermi energies. The calculations were performed 
on a 50 X 50 lattice with the magnetic field and disorder turned 
off. The table also shows the exact values of Gu obtained with 
the classical Kubo formula and the relative errors ocurring in 
the first calculation. 



torus until the computation becomes virtually exact. At 
the same time, we can still evaluate the conductivity via 
Eq. [72] A comparison between the two results will pro- 
videthe test. 

Table [l] lists an as function of Fermi energy, with the 
noncommutative Kubo formula of Eq.|72]evaluated on a 
50 X 50 lattice. The table also lists the virtually exact an 
values obtained with the traditional Kubo formula in fc- 
space. As the table shows, all the errors are less than one 
percent, except in one case where we see a 2% error (due 
to a spike in the density of states at that energy). Based 
on this table we expect well converged results when the 
calculations are performed on 50 x 50. When disorder 
is present, the convergence will be double-checked by 
comparing the outputs from calculations performed on 
lattices of increasing sizes, typically 30 x 30, 40 x 40 and 
50 X 50. As we shall see, the convergence of the calcula- 
tions is also confirmed by those tests. 



C. The diagonal conductivity as function of disorder and 
Fermi level 

Fig. [T] reports the calculated an for the topological 
case M = 6 (the off-diagonal components of a are zero 
because of time reversal symmetry). The diagonal con- 
ductivity is plotted as function of Fermi energy at various 
disorder strengths. The first three columns show the re- 
sults obtained by evaluating Eq. [72]on 30 X 30, 40 x 40 
and 50 x 50 lattices. The temperature and the relaxation 
time were fixed at: fcT = l/xrei = 0.01. The calculations 
were repeated for 10 random disorder configurations 
and the results are all shown in the first three columns of 
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FIG. 1. The diagonal conductivity for the topological case (M = 6) as function of Fermi level, lattice size and disorder strength. For 
each Fermi level, the noncommutative Kubo formula was evaluated for 10 random disorder configurations, with kT = 1/ z^ei = 0.01. 
The last column shows the averages over the disorder configurations. The averages corresponding to the three lattice sizes overlap 
each other almost perfectly indicating a good convergence of the calculations with the lattice size. 



Fig. [T] w^ithout averaging (hence the fuzziness displayed 
in those graphs). The (vertical) spread of an due to the 
changing of the disorder configuration tells us about the 
degree of self -aver aging in these calculations. Clearly 
the spread decreases, and as such, the self-averaging is 
improved as the lattice size is increased. But the most im- 
portant fact about these data, are the averages of on over 
the disorder configurations shown in the 4-th column of 
Fig. [T] There are minor or no visible variations between 
the averages obtained on different lattice sizes, a fact 
that indicates a very well converged simulation. With 
that, we fulfilled the main goal of this work, namely, to 
demonstrate the efficiency and effectiveness of the non- 
commutative Kubo formula. 

We now discuss the data. In Fig.[lj one can see that, if 
the Fermi level is located in certain energy regions, an 
rapidly decreases as the disorder strength is increased. 
But one can also see well defined energy regions where 
On saturates and stays at a certain appreciable value 
even at large disorder strengths. As we shall see from 



the temperature-dependence analysis, the latter spectral 
regions are metallic in character, while the former ones 
are insulating. 

Fig. [ij also displays a pattern that is very specific to 
topological models. The energy regions where on is 
maximum are seen to drift towards each other and ul- 
timately to merge together as the disorder strength is 
increased. This phenomenon is called the levitation 
and annihilation of the conducting states. Its origin is 
connected to the fact that the conducting states carry a 
non-trivial topological Z2 invariant. Since this invariant 
is robust against strong disorder, the conducting states 
cannot suddenly disappear and, instead, neighbouring 
conducting states levitate towards each other and anni- 
hilate their Z2 invariants. From there on, the states can 
localize, as it inherently happens in any system if the 
disorder strength is large enough.^^ This phenomenon 
was widely discussed in the context of Integer Quantum 
Hall effect,^^*^^^and it was first observed in the context 
of Quantum spin-Hall Insulators in Ref..68iand in Chern 
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FIG. 2. The diagonal conductivity for the topological case 
(M = 6) as function of Fermi energy, disorder strength and 
temperature. An average over 10 disorder configurations was 
used. Each panel displays three curves, corresponding to 
kJ = 1/Trei = 0.01 (black), kT = l/z^ei = 0.025 (gray) and 
kT = l/Xrei = 0.05 (light gray). The shaded regions indicate 
the Fermi energies where On increases when the temperature 
is reduced, i.e. where the model displays a metallic behavior. 



insulators in Ref . [ 



D. Temperature dependence and the phase diagram 

The relaxation time Trei behaves as Trei ~ l/(fcT)" (a > 
0) in the limit of low temperatures, where the exponent a 
depends on the dominant dissipation mechanism (a = 5 
if the dissipation is through phonons). As such, the 
temperature dependence of an comes from the Fermi- 
Dirac statistics and from Trei- The behavior of on as 
function of T in the asymptotic limit T ^ can provide 
an accurate picture of the nature of the energy spectrum 
and the transport characteristics of the system. 



The spectral, fractal and diffusion exponents and their 
inter-relations and relevance to transport in aperiodic 
solids were introduced in Ref. 40 In particular, the 
diffusion exponent jSdiff at the Fermi level was defined 
through the mean square of the displacement operator 
in the asymptotic limit of large time intervals: 

dx^{t) = - f dt' r({x{t') - x)^nE,) f^^^'. (79) 
^ Jo 

Here x is the position operator, x{t) is its time evolution 
and 7Tf^ is the projector onto the energy spectrum below 
Ef . The following behavior of on at low temperatures 
was proved in Ref. I40l 



On 



J«(l-2i6di 



(80) 



This relation is the Drude formula for aperiodic solids. 

The diffusion exponent takes values between and 1. 
Its precise value can give information about the nature 
of energy spectrum. In general, the energy spectrum can 
be absolute-continuous, singular-continuous (fractal) or 
pure point (localized). The transport is called ballistic 
if i^diff = 1 and that requires the spectrum around Ef to 
be absolute continuous. The reverse statement is true in 
one dimension, but in higher dimensions there are mod- 
els with absolute continuous spectrum but with non- 
ballistic (diffusive) transport.^^ The transport is called 
diffusive if < jSdiff < 1 and occurs when the energy 
spectrum around Ef is singular-continuous. Absence 
of diffusion jSdiff = occurs when the spectrum is pure 
point (localized). A detailed discussion of jSdiff for differ- 
ent quantum models, together with the physical implica- 
tions, was given in Refs. 70 and 52 We want to point out 
that the exponent in Eq.[8Q|was experimentally evaluated 
in Ref. 25 for a film of Bi2Se3. We believe this exponent 
can be extracted from many other experimental data. 

While we do have the technology to compute the dif- 
fusion exponent for our topological model, here we ad- 
dress a coarser problem, namely, we determine when the 
system behaves as a conductor: an y oo as T \ 0, or as 
an insulator: an \ as T \ . This is something we 
can do at the grand-scale of Fig.llJ which then will allow 
us to draw the phase diagram orthe system. A study on 
the jSdiff itself will be reported elsewhere. 

For this we repeated the calculations from Fig. [T] for: 
kT = 1/Trei = 0.025 and kT = l/xrei = 0.05. We thus 
assume an exponent a = 1 but note that the conducting 
or insulating character is independent of this exponent. 
The theory predicts a much faster convergence for these 
cases towards the thermodynamic limit, so we restricted 
these calculations to only a 40 x 40 lattice. The results 
are reported in Fig. |2] According to the previous dis- 
cussion, in these plots we should see conducting energy 
regions where the values of an increase as T is lowered 
and insulating energy regions where an decreases as T 
is lowered. This is indeed consistently observed in Fig.|2] 
for all three curves corresponding to different tempera- 
tures. For example, all three curves intersect each other 
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FIG. 3. The diagonal conductivity for the trivial case (M = 9) as function of Fermi level, lattice size and disorder strength. For each 
Fermi level, the noncommutative Kubo formula was evaluated for 10 random disorder configurations, with /cT = l/Xrei = 0.01. 
The last column shows the averages over the disorder configurations. The averages corresponding to the three lattice sizes overlap 
each other almost perfectly indicating a good convergence of the calculations with the lattice size (maybe with the exception of 
the last panel). 



at more or less the same point, so the transition point 
between the conducting and insulating regions can be 
identified. As expected, the conducting energy regions 
occur in the middle of the bands while the insulating 
regions occur near the edges. In agreement with the lev- 
itation and annihilation picture, the conducting regions 
are seen to drift towards each other until they merge, at 
about W = 7, at which point they rapidly diminish as W 
is being further increased. 



E. Disorder induced conducting states 

The model used in our simulations can enter the 
Quantum spin-Hall topological phase upon increasing 
the disord er, eve n if we start from the trivial phase 
at W = Qj 45 | 62 | 7iE | jg (j^g |-Q ^ strong disorder- 
induced deformation of the Quantum spin-Hall phase 
boundary.^^ We should mention that this phenomenon 



can disappear if other types of disorder are used.^^ The 
phase diagram of the model with onsite disorder was 
computed by various methods such as by mapping 
the conductance of the edge states,^^^ the Lyapunov 
exponent,^^ or the spin-Chern number.HS 

With the Rashba interaction turned on, the model be- 
longs to the symplectic class so the topological and trivial 
phases are separated by a metallic phase.^^^^ As such, 
we can use our transport simulations to see if indeed 
there are metallic states induced at large disorder, when 
such states are absent at weak disorder. In principle, 
we can use the transport simulations to map the whole 
metallic phase between the trivial and the topological 
insulating phases, in the 3-dimensional space of £f, W 
and M, and that will be the fourth way to compute the 
phase diagram of the model. The fifth way will be to 
use the level statistics analysis as it was do ne for ot her 
models of disordered topological insulators .ESSEZIZSI 

However, in this section we will compute just another 
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slice of the phase diagram, corresponding to M = 9 
which is in the trivial part of phase diagram at W = 0. We 
have repeated the simulations reported in Fig.[T](with M 
set at the new value) and the results are reportea in Fig.lS] 
Here, we see again a reduction in the spread due to the 
disorder of on as the size of the lattice is increased, and 
a good overlap of the disorder-averages for the different 
lattice sizes. While for the most part of the spectrum the 
conductivity rapidly decay with the increase of W, one 
can see in the middle of the spectrum and starting from 
W = 7 a sudden increase of an until it reaches the same 
value as in Fig.jl] To demonstrate that the system enters 
a metallic phase, we have repeated the simulations from 
ing M the new value. The results are reported in 
[ they indeed confirm that, starting from W = 7, 
there is a region where an increases as the temperature 
is lowered, hence those states have metallic character. 

We will like to point out that our conclusions based on 
the transport calculations are in quantitative agreement 
with the phase diagrams computed in Refs. [62] and |45l 
For example, our data in Fig. [l] show that at W = 7 
the topological phase is completely gone, and only the 
metallic and the trivial phases are still present. This is 
precisely what was predicted in these two mentioned 
studies. Furthermore, it was predicted that if M = 9, 
then the model will enter the metallic phase at about 
W = 7 which is exactly what we observe in our transport 
simulations. 



VII. COMPARISON WITH THE LEVEL STATISTICS 

A. The topological case 

We performed a level statistics analysis in precisely 
the same manner as in our previous studies. ^'^ The 
exact procedure for the level statistics analysis was de- 
scribed in detail in these publications. Fig. |5] reports the 
variance of the distribution of the level spacings collected 
from small energy windows centered at various energies 
while the random potential was updated 500 times. The 
mass term was fixed at M = 6. It is a well established fact 
that/^ if the localization length of the system is smaller 
than the simulation box, then the level spacings follow a 
Poisson distribution which has variance equal to 1, and 
if the localization length is comparable or larger than the 
simulation box, then the level spacings follow the Wigner 
surmise for symplectic Gaussian ensembles which has a 
variance of 0.104. As such, the Fig.|5]allows us to identify 
the spectral regions with very large or infinite localiza- 
tion lengths, which coincide with the energy intervals 
where the variance is close to expected value of 0.104. 
The level statistics analysis was performed on a 40x40 
lattice. 

The energy regions harboring the extended states have 
been shaded in Fig.|5]for a better visualization. The phe- 
nomenon of levitation and annihilation is clearly dis- 
played there and the emerging phase diagram is in good 
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FIG. 4. The diagonal conductivity for the trivial case (M = 9) as 
function of Fermi energy, disorder strength and temperature. 
An average over 10 disorder configurations was used. Each 
panel displays three curves, corresponding to kT = l/Xrei = 
0.01 (black), kT = l/T^ei = 0.025 (gray) and kT = l/Xrei = 0.05 
(light gray). The shaded regions indicate the Fermi ener- 
gies where Ou increases when the temperature is reduced, i.e. 
where the model displays a metallic behavior. 



quantitative agreement with Fig. [2] Strictly speaking, 
the metallic region identified in Fig. |2] is strictly con- 
tained (thus not equal) inside the region of extended 
states identified in Fig. |5] That is because the metallic 
phase in Fig. |2] contains all the states with jSdiff > 0.5 
while the phase of extended states in Fig. |5] contains all 
the states with jSdiff > 0. 



B. The disorder induced conducting states 

We have repeated the level statistics analysis f or M = 9 
and the results are reported in Fig. [6] We have again 
shaded the energy regions that harbor extended states 
and, sure enough, we observe the emergence of the ex- 
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FIG. 5. The variance of the ensembles of level spacings for the 
topological case (M = 6), recorded at various energies as func- 
tion of disorder strength. A total of 500 disorder configurations 
were used in these simulations and, for each disorder configu- 
ration and energy E, 13 level spacings were collected from the 
immediate vicinity of £. As such, the ensembles contain 6,500 
level spacings. The size of the lattice for these simulations was 
40 X 40. 
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FIG. 6. Same as Fig. |5] for the trivial case (M = 9). 



states and the extended states in Fig.|5]behave as the dis- 
order is increased. In the former case, the levitation and 
annihilation is absent and the energy domains harboring 
the extended states simply contract and vanish. 



tended states that were previously seen in our transport 
simulations reported in Fig. [4] Besides these interesting 
conducting states, one can also see the extended states 
that originate directly from the bands of the clean system. 
We should probably mention that in a 2-dimensional 
symplectic model, unlike the unitary Anderson model, 
the extended states can survive moderate disorder. But 
there is a distinct difference in the way these extended 



VIII. TRANSPORT SIMULATIONS FOR DISORDERED 
QUANTUM SPIN-HALL INSULATORS UNDER 
MAGNETIC FIELDS 

The time-reversal symmetry is essential for the topo- 
logical properties of the Quantum spin-Hall Insulators. 
Nevertheless, the transport experiments in the presence 
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Energy 

FIG. 7. Evolution of the density of states (DOS) for the topo- 
logical case M = 6 (left) and trivial case M = 9 (right) with 
the variation of the magnetic field. The lattice size for this 
simulations was 100x100 and W = 0. 



of a magnetic field, which breaks the time-reversal sym- 
metry, are extremely useful for understanding the mi- 
croscopic electronic structure of the samples. In a typ- 
ical experiment, the transport coefficients are mapped 
as functions of the magnetic field strength and of a gate 
potential. The Hall resistivity measurements at weak 
magnetic fields provide an accurate map of the electron 
density as function of the gate potential, and typically 
this function has a simple linear shape. As such, the de- 
pendance of the transport coefficients on the magnetic 
field, electron density and temperature can be accurately 
determined. 

A typical transport measurement on HgTe Quantum 
spin-Hall insulators or thin films of 3-dimensional topo- 
logical insulators shows the following qualitative fea- 
tures: Initially, the Hall resistivity increases (decrease) 
linearly for n-type (p-type) carrier concentration as the 
magnetic field is turned on. Hall plateaus start to ap- 
pear at some point and they become fairly wide and well 
defined at larger field strengths, typically of a few Tes- 
las. The diagonal resistivity starts flat as the magnetic 
field is strengthened, but at some point it starts to de- 
velop well defined oscillations of increasing amplitude, 
the Shubnikov-de Haas oscillations. At strong magnetic 
fields, where the Hall plateaus are wide and well de- 
fined, the direct resistivity becomes zero inside the Hall 
plateaus and displays sharp peaks at the edges of the 
Hall plateaus. The Shubnikov-de Haas oscillations can 
be analyzed using semi-classical theoretical approaches 
and information about the geometry of the Fermi sur- 
face can be extracted. This type of analysis has been 
intensely used in the search of the 2-dimensional com- 
ponent of the Fermi surface of 3-dimensional topological 
insulators, which can be related to the topological surface 
states. Besides the Shubnikov-de Haas oscillations, the 
widths of the Hall plateaus could be used to assess the 
degree of disorder present in the samples. Our hope is 
that the present technique will enable us to reverse engi- 
neer the experimental results on the resistivity tensor to 
a point where we can pin-point a particular microscopic 
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FIG. 8. Upper panel: The diagonal and Hall resistivities for the 
topological case (M = 6), computed at a fixed magnetic flux 
(j) = 0.08 X <po while the Fermi energy was swept over the entire 
energy spectrum. The lattice size in these simulations was 
50x50, W = and /cT = l/T^ei = 0.01. Lower panel: The data 
from the upper panel is replotted as function of the electron 
density. 



model for a given sample. 

The transport simulations in the presence magnetic 
fields present the same technical difficulties as before 
since the effect of the magnetic fields is introduced via 
the Peierls substitution as discussed. However, to re- 
solve all the expected features, we will need to increase 
the simulation box. In this section we only want to 
demonstrate that we can indeed perform such simula- 
tions using the non-commutative formalism, more pre- 
cisely, that we can resolve Hall plateaus and Shubnikov 
de Haas oscillations when the resistivity tensor is plotted 
as function of electron density. These features are easily 
reproduced when the data is plotted as function of the 
Fermi energy, but as we shall see, they disappear when 
plotted as function of electron density (which is how the 
data is plotted in the experiments), unless fairly strong 
disorder is present. As such, resolving the Hall plateaus 
in a transport simulation is a highly non-trivial problem, 
which even for the simple Integer Quantum Hall Effect 
has not been completely solved. 

We now start the discussion of our simulation results. 
The energy spectrum of periodic systems under mag- 
netic fields have a Hofstadter fractal structure.'''' Fig.jT^ 
shows the density of states (DOS) as function of energy 
and magnetic field, when the model is in the topological 
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FIG. 9. Same as Fig.[7|but for W = 3. 



and trivial phases, with the disorder turned off. As it 
was already pointed out in Ref. 13. there is an interest- 
ing qualitative difference in the behavior of the energy 
spectrum for the topological and trivial cases. What one 
basically sees in Fig.[7|is two Hofstadter butterflies, one 
below and one above the spectral gap, that are attracted 
to each other in the topological case, and are repelled 
by each other in the trivial case. The behavior of DOS 
with the magnetic field reminds us of the levitation and 
annihilation phenomenon discussed in the previous sec- 
tions, but the cause of this behavior is not understood 
yet. Given that the insulating gap of most topological 
materials is small, this effect must be definitely taken 
into account when interpreting the experimental data. 

Now the Hofstadter spectrum is complex and has a 
fractal nature, but still at low magnetic fields one can 
distinguish clear thin Landau bands in Fig.|7| Of course, 
if we zoom in on these bands, we will see the whole 
structure repeating itself over and over again. But at 
the scale of Fig. [Zj when the Fermi level is located in 
between these Landau bands, the Hall conductivity or 
resistivity should be quantized and the direct conduc- 
tivity or resistivity should display a local minimum (the 
simulations are at finite temperature so we don't expect 
a strictly null direct conductivity). The Hall conductiv- 
ity or resistivity should jump to the next Hall plateau 
when the Fermi level is brushed over a Landau band, 
while the direct conductivity or resistivity should have 
a spike. The Landau bands are clearly contoured at the 
bottom or at the top of the energy spectrum but, at the 
scale of Fig.[7| they are not well resolved at the edges of 
the insulating gap where one will be most interested in. 
As such, our transport simulations in which the Fermi 
level is varied over the entire energy spectrum and are 
done at the same scale as that of Fig. U\ are expected to 
show quantized Hall plateaus only at the bottom and at 
the top of the energy spectrum. More refined simula- 
tions which hopefully could resolve the energy region 
near the gap edges will be presented elsewhere. 

Fig.jSjreports the simulated direct and Hall resistivities 
for a magnetic flux per unit cell oi(p = 0.08 (po. The sim- 
ulations were performed on a 50 x 50 lattice. In Fig.|8|a) 
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FIG. 10. Same as Fig.|8]but for W = 3. 



the results are plotted as function of Fermi energy, while 
in Fig.[8|b) the results are plotted as function of electron 
density. In the first case, the Hall plateaus appear very 
clearly contoured and the spikes in the direct resistivity 
can be seen whenever the sampled energies fell between 
two Hall plateaus. However, when the data is plotted 
as function of electron density, the Hall plateaus are re- 
duced to a single point because the electron density does 
not vary when the Fermi level sweeps over the clean 
spectral gaps. As a result, the Hall plateaus disappear. 
When the disorder is turned on, the Hofstadter spectrum 
is smoothed out by the localized states that are pulled 
out the Landau bands, as one can see in Fig. [9] where we 
plot the DOS for W = 3. The results of the transport sim- 
ulations for W = 3 are reported in Fig.jlO] The magnetic 
flux per unit cell and the lattice size were fixed at the 
same values as for Fig.|8] The plots in Fig.[9]are obtained 
with a single disorder configuration. To read the data 
in Fig. [9j it is best to look first at the direct conductivity 
because it displays four clear dips. Examining the Hall 
resistivity plots, we see that they display a Hall plateau 
at each of these dips. The Hall plateaus are not perfectly 
quantized but their widths are clearly identifiable and 
that is enough for comparisons with the experiments. 
Very importantly, the widths of the Hall plateaus remain 
finite when the data is plotted as function of the electron 
density. 
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IX. CONCLUSIONS 



The main purpose of the present work was to intro- 
duce the reader to the noncommutative theory of charge 
transport and to present an efficient and fast converg- 
ing numerical implementation of the noncommutative 
Kubo formula. As an interesting application, we chose a 
model of a topological insulator where we were able to 
map the diagonal conductivity as function of Fermi level, 
disorder strength and temperature. This enabled us to 
demonstrate that the topological phase is surrounded 
by a diffusive metallic phase where the transport co- 
efficients increase as the temperature is lowered. This 
means that disorder alone can drive the system from a 
topological insulator into a diffusive metal. We postulate 



that this is what is actually happening in the current ex- 
perimental observations where all the samples showed 
metallic bulk so far. When we introduce a perpendicu- 
lar magnetic field, we found an intriguing behavior of 
the energy spectrum in that the valence and conduc- 
tion bands move towards each other as the strength of 
the field is increased, until they touch and merge. In 
our transport simulations at fixed magnetic flux and in 
the presence of disorder, we were able to resolve a few 
Hall plateaus even when the conductivity was plotted 
as function of electron density. 
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